# Examine how the quality of liberal democracy conditions the effects of populist involvement on post-election 
# emotions. See Table N.1 and Figure N.6, Appendix N.




#############
# load data #
#############

library(brms)
library(ggplot2)
library(ggpubr)
library(stargazer)

load("PSRM Replication Files/PresidentialElectionsPosts.RData")

colnames(president_df)

all(c("populist_present", "v2x_libdem") %in% colnames(president_df)) # TRUE

length(unique(president_df$country)) # 26
length(unique(president_df$electionname)) # 29
length(unique(president_df$party_country)) # 52

# focus on 15 days before and after the election
day_range <- 15

sub <- president_df[which(president_df$daysinceelection >= -1*day_range
                          & president_df$daysinceelection <= day_range),]

# winners and losers
winners <- sub[which(sub$win == 1),]
losers <- sub[which(sub$win != 1),]

length(unique(winners$electionname[which(!is.na(winners$polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$economic_polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$social_polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$populist_present))])) # 27
length(unique(winners$electionname[which(!is.na(winners$populist_dummy))])) # 27

length(unique(losers$electionname[which(!is.na(losers$polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$economic_polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$social_polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$populist_present))])) # 27
length(unique(losers$electionname[which(!is.na(losers$populist_dummy))])) # 27




##################
# winners + love #
##################

winner_love_qd <- brm(loveprop ~ post*populist_present*v2x_libdem
                      + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                      + Xl + Xr 
                      + (1 + Xl + Xr | party_election),
                      data=winners,
                      warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(winner_love_qd)[[14]][,5])
#plot(winner_love_qd)
#summary(winner_love_qd)




###################
# winners + angry #
###################

winner_angry_qd <- brm(angryprop ~ post*populist_present*v2x_libdem
                       + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                       + Xl + Xr 
                       + (1 + Xl + Xr | party_election),
                       data=winners,
                       warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(winner_angry_qd)[[14]][,5])
#plot(winner_angry_qd)
#summary(winner_angry_qd)




#################
# losers + love #
#################

loser_love_qd <- brm(loveprop ~ post*populist_present*v2x_libdem
                     + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                     + Xl + Xr 
                     + (1 + Xl + Xr | party_election),
                     data=losers,
                     warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(loser_love_qd)[[14]][,5])
#plot(loser_love_qd)
#summary(loser_love_qd)




##################
# losers + angry #
##################

loser_angry_qd <- brm(angryprop ~ post*populist_present*v2x_libdem
                      + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                      + Xl + Xr 
                      + (1 + Xl + Xr | party_election),
                      data=losers,
                      warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(loser_angry_qd)[[14]][,5])
#plot(loser_angry_qd)
#summary(loser_angry_qd)

#save(winner_love_qd,
#     winner_angry_qd,
#     loser_love_qd,
#     loser_angry_qd,
#     file="Stan_PopulismQuaityDemocracy.RData")




##################
# result summary #
##################

load("PSRM Replication Files/Stan_PopulismQuaityDemocracy.RData")

round(posterior_interval(winner_love_qd, "post", 0.95), 2)
round(posterior_interval(winner_angry_qd, "post", 0.95), 2)
round(posterior_interval(winner_love_qd, "post", 0.95), 2)
round(posterior_interval(loser_angry_qd, "post", 0.9), 2)

# table
f1 <- lm(loveprop ~ post*populist_present*v2x_libdem
         + incumbentparty + concurrent + runoff + semipresidential + enc_v1
         + Xl + Xr,
         data=winners)

f2 <- lm(angryprop ~ post*populist_present*v2x_libdem
         + incumbentparty + concurrent + runoff + semipresidential + enc_v1
         + Xl + Xr,
         data=winners)

f3 <- lm(loveprop ~ post*populist_present*v2x_libdem
         + incumbentparty + concurrent + runoff + semipresidential + enc_v1
         + Xl + Xr,
         data=losers)

f4 <- lm(angryprop ~ post*populist_present*v2x_libdem
         + incumbentparty + concurrent + runoff + semipresidential + enc_v1
         + Xl + Xr,
         data=losers)

# table N.1
stargazer(f1, f2, f3, f4,
          coef=list(fixef(winner_love_qd)[,1], fixef(winner_angry_qd)[,1], 
                    fixef(loser_love_qd)[,1], fixef(loser_angry_qd)[,1]),
          #se=list(fixef(winner_love_qd)[,2], fixef(winner_angry_qd)[,2], 
          #        fixef(loser_love_qd)[,2], fixef(loser_angry_qd)[,2]),
          ci.custom=list(fixef(winner_love_qd)[,3:4], fixef(winner_angry_qd)[,3:4], 
                         fixef(loser_love_qd)[,3:4], fixef(loser_angry_qd)[,3:4]),
          omit="Constant",
          star.cutoffs=NA, 
          digits=2,
          no.space=TRUE,
          dep.var.caption="",
          omit.stat = c("f", "adj.rsq", "rsq", "ser"),
          dep.var.labels=c("Love", "Angry", "Love", "Angry"),
          covariate.labels=c("Post Election",
                             "Populist Involvement", "Quality of Liberal Democracy",
                             "Incumbent Party", "Concurrent Election",
                             "Runoff", "Semi-Presidential",
                             "Effective Number of Candidates",
                             "Pre-Election Trend (Group Mean)", 
                             "Post-Election Trend (Group Mean)",
                             "Post Election times Populist Involvement",
                             "Post Election times Quality of Liberal Democracy",
                             "Populist Involvement times Quality of Liberal Democracy",
                             "Post Election times Populist Involvement times Quality of Liberal Democracy"))

# standard deviation of random effects
getSigmas <- function(model){
  sest <- round(summary(model)[[17]]$party_election[1:3,1], 2)
  ssd <- paste0("[",
                sapply(1:3, function(p){paste(round(summary(model)[[17]]$party_election[p,3:4], 2), 
                                              collapse = ", ")}), "]")
  
  return(c(sest[1], ssd[1], sest[2], ssd[2], sest[3], ssd[3]))
}

sigmas1 <- data.frame(label=c("hatsigmatextIntercept", "",
                              "hatsigmatextPre-Election Trend", "",
                              "hatsigmatextPost-Election Trend", ""),
                      s1=getSigmas(winner_love_qd),
                      s2=getSigmas(winner_angry_qd),
                      s3=getSigmas(loser_love_qd),
                      s4=getSigmas(loser_angry_qd))
stargazer(t(t(sigmas1)))

c(length(unique(winner_love_qd$data$party_election)),
  length(unique(winner_angry_qd$data$party_election)),
  length(unique(loser_love_qd$data$party_election)),
  length(unique(loser_angry_qd$data$party_election)))


# plot marginal effects
getME <- function(model, populist_present, v2x_libdem, label){
  posteriors <- posterior_samples(model, "b_post")
  
  est <- data.frame(t(sapply(1:length(v2x_libdem), function(x){
    post_est <- unlist(posteriors[,1] + posteriors[,2]*populist_present + posteriors[,3]*v2x_libdem[x] + posteriors[,4]*populist_present*v2x_libdem[x])
    return(c(populist_present,
             v2x_libdem[x], 
             mean(post_est),
             quantile(post_est, 0.025),
             quantile(post_est, 0.975)))
  })))
  
  colnames(est) <- c("populism", "libdem", "est", "lwr", "upr")
  
  est$label <- label
  
  return(est)
}

range(president_df$v2x_libdem)
libdem_sim <- seq(0.3, 0.85, length=100)

metab <- rbind(getME(winner_love_qd, 0, libdem_sim, "(a) Winner + Love"),
               getME(winner_love_qd, 1, libdem_sim, "(a) Winner + Love"),
               getME(winner_angry_qd, 0, libdem_sim, "(b) Winner + Angry"),
               getME(winner_angry_qd, 1, libdem_sim, "(b) Winner + Angry"),
               getME(loser_love_qd, 0, libdem_sim, "(c) Loser + Love"),
               getME(loser_love_qd, 1, libdem_sim, "(c) Loser + Love"),
               getME(loser_angry_qd, 0, libdem_sim, "(d) Loser + Angry"),
               getME(loser_angry_qd, 1, libdem_sim, "(d) Loser + Angry"))

metab$populism <- ifelse(metab$populism == 1, "Yes", "No")

metab$label <- factor(metab$label, levels=unique(metab$label))

# figure N.6
ggplot(data = metab, 
       aes(x = libdem, y = est, color = populism, fill = populism)) + 
  geom_ribbon(data = metab, aes(ymin = lwr, ymax = upr), alpha = 0.3, linetype = 0) +
  geom_hline(yintercept = 0, lwd = 0.5, lty = 2) +
  geom_line(lwd=1.5) +
  scale_color_manual(name = "Populist Involvement", values = c("#0a8dff", "#ff640a")) +
  scale_fill_manual(name = "Populist Involvement", values = c("#0a8dff", "#ff640a")) +
  facet_wrap(. ~ label, ncol = 2) +
  xlab("Quality of Liberal Democracy") +
  ylab("Marginal Effect of Post Election") +
  theme_bw() +
  guides(color=guide_legend(override.aes=list(fill=NA))) +
  theme(legend.position = "bottom",
        strip.text.x = element_text(size = 14),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 16),
        axis.title.y = element_text(size = 16))

#ggsave("marginaleffect_qualitydemocracy.pdf", width=10, height=10)
